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Abstract 

We investigate the volume growth of ionized regions around UV photon sources 
with the WENO algorithm, which is an effective solver of photon kinetics in the 
phase space described by the radiative transfer equation. We show that the volume 
growth rate, either of isolated ionized regions or of clustered regions in merging, 
generally consists of three phases: fast or relativistic growth phase at the early 
stage, slow growth phase at the later stage, and a transition phase between the 
fast and slow phases. The growth rate can be characterized by a time scale tc of 
the transition phase, which is approximately proportional to E^^'^^ E being the 
intensity of the ionizing source. The larger the time scale tc, the longer the photons 
to postpone their contribution to the ionization. For strong sources, like E > 10^^ erg 
s~^, tc can be as large as a few Myrs, which could even be larger than the lifetime 
of the sources. Consequently, most photons from these sources contribute to the 
reionization only when these sources already ceased. We also show that the volume 
growth of ionized regions around clustered sources with intensity Ei [i = 1,2, . . .) 
would have the same behavior as a single source with intensity E = Ei , if all the 
distances between nearest neighbor sources i and j are smaller than c(t* + ti), 
being the time scale tc of source i. Therefore, a tightly clustered UV photon sources 
would lead to a slow growth of ionized volume. This effect would be important for 
studying the redshift-dependence of 21cm signals from the reionization epoch. We 
also developed, in this paper, the method of using WENO scheme to solve radiative 
transfer equation beyond one physical dimension. This method can be used for high 
dimensional problems in general as well. 
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1 Introduction 



In the early stage of reionization of the universe, radiation from the first 
generation of stars ionizes neutron hydrogen and hehum to produce ionized 
bubbles around these stars. The subsequent growing, overlapping and merging 
of these isolated ionized patches lead to a full reionization of the universe. The 
evolution of the reionization depends on the birth rate of the first stars and the 
formation of ionized regions around these sources. The growth of ionized HII 
volume is directly related to the formation of 21 cm emission and absorption 
regions at the reionization epoch (e.g. Cen 2006; Alvarez et al. 2006; Chuzoy 
et al. 2006; Liu et al. 2007). The planning and ongoing projects of detecting 
redshifted 21 cm signals are trying to reconstruct the redshift evolution of 
the reionized volume. Therefore, a detailed study on the merging of ionized 
regions is necessary. 

The growth of the ionized HII volume, V, is usually described by a rate equa- 
tion (e.g. Shapiro & Giroux 1987; Madau, Haardt & Rees 1999; Mellema et 
al. 2006): 

n{t)^ = TV - / n-'{t)aBC{t)dV, (1) 



where n = 1.88 x 10"^(r2fe/i^/0.022)(l + z)^ cm~^ is the mean number density 
of hydrogen at redshift z, N is the emission rate of ionizing photons, as is the 
recombination coefficient, and C{t) is the volume- averaged clumping factor of 
HII. 

From equation (1), V{t) is linearly dependent on N: V{t) oc N. That is, the 
growth of the ionized volume is proportional to the emission rate of ionizing 
photons N, regardless of whether the ionizing photons are produced from one 
source with emission rate N or from m sources with emission rate N/m. The 
linear relation between V(t) and N has been used in the simulations of the 
reionization (e.g. Ciardi et al. 2003; Iliev et al. 2006). That is, the effects of 
the photon kinetics in the phase space are completely ignored. 

A basic assumption of equation (1) is that photons emitted from the sources 
will immediately join the action with the atom, regardless of the photon prop- 
agation between the source and the atom. This assumption is reasonable if 
the time scale of the photon kinetics is much smaller than that of the prob- 
lem considered. Unfortunately, this is not always correct for problems at the 
epoch of reionization (Shapiro et al. 2006; Qiu et al. 2007). In this paper, 
we will show that the finite speed of light will lead to a substantial change 
of the growth rate of the ionized volume. The ionized volume growth of one 
source with emission rate N can be significantly different from that of multiple 
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sources with the total emission rate equal to N. The growth rate of the ionized 
volume depends not only on the total emission rate of ionizing photons, N, 
but also on the distribution and clustering of the sources. 

Many numerical solvers for the radiative transfer equation have been proposed 
(Razoumov & Scott 1999; Abel et al. 1999; Ciardi et al. 2001, Gnedin & Abel 

2001, Sokasian et al. 2001, Nakamoto et al. 2001; Razoumov et al. 2002, Cen 

2002, Maselh et al. 2003, Shapiro et al., 2004; Rijkhorst et al. 2006; Mellema 
ct al. 2006; Susa 2006, Whalen & Norman 2006). We use the WENO scheme 
to be the solver for the photon kinetics in the phase space. The WENO al- 
gorithm has been proved to have high order accuracy and good convergence 
in capturing discontinuities and complicated structures in fluid as well as to 
be significantly superior over piecewise smooth solutions containing discon- 
tinuities (Shu 2003). We have showed that the WENO algorithm is effective 
for solving radiative transfer problem in one-dimensional physical space and 
one-dimensional frequency space (Qiu et al. 2006, 2007). It revealed that the 
time-dependent solution of the radiative transfer equations is essential for the 
formation and evolution of the ionized and heated regions around UV ionizing 
sources. In this paper, we will develop the WENO algorithm of the radiative 
transfer equations beyond one-dimension. 

The paper is organized as follows. Section 2 presents the source-intensity de- 
pendence of the growth rate of ionized regions around isolated point sources. 
Section 3 studies the merging of two ionized regions and its effect on the growth 
of the ionized volume. Discussions and conclusions arc given in Section 4. The 
details of the WENO numerical scheme are listed in Appendix. 



2 Isolated point source 



The patchy structures of the HII region in the early universe is very com- 
plex. However, in the early stage, many HII regions are isolated and even 
spherical around UV photon sources. Later, these spherical regions merge and 
yield complicated structures. As a preparation, we summarize in this section, 
the features of the growth of the isolated ionized regions, calculated with the 
WENO algorithm (Qiu et al. 2006, 2007, Liu et al. 2007). To make the pa- 
per self-contained, the corresponding equations and parameters are given in 
Appendix A. 
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2. 1 The growth of the ionized volume 



In deriving eq.(l), it is assumed that the ionization in region V is perfect, i.e., 
the fraction of the neutral hydrogen, /hi = nni/nn, is zero within the region, 
and 1 outside. Actually, the ionization cannot be complete, and the ionization 
front (I-front) cannot be defined by a sharp boundary dividing the completely 
ionized and completely neutral regions. We will define the ionized region to 
be the place in which /hi < 90%. 

Consider a point source emitting photons of energy E{i/)di' per unit time 
within the frequency ranges from i/ to I' + di'. The energy spectrum of photons 
is assumed to be a power law, i?(z/) = Eo{i>o/i')°' with a = 2, uq is the 
ionization energy. Assuming the hydrogen gas around the source is uniform, 
the volume growth of the ionized region, V{t), at redshift 1 + z = 10, is 
shown in Figure 1, in which the intensities of UV photons are taken to be 
E = /~ E{i^)di^ = 5.8 X 10=^^ 10^^ 10^3 and 10^^ erg s-^ or TV = 1.34 x 10™, 
10^^ 10^^ and 10™ s'K 

In Figure 1, we use Myrs and Mpc to be the units of time t and length r. In 
this paper, we also sometimes use dimensionless time and length defined as 
t' = cnaot and r' = na^r^ where ctq is the ionization cross section. Therefore, 
t = 0.89(1 + zyH' Myrs and r = 0.27(1 + z)~^r' Mpc. t' and r' actually are in 
the units of mean free fiight time and mean free path of ionizing photons. The 
dimensionless variables are convenient for numerical works (see Appendix C). 

Figure 1 shows a common feature for all sources that the growth of the ionized 
volume undergoes three phases: when t is small, the growth is very fast, when 
t is large the growth is very slow, and a transition phase between the fast and 
slow phases. The details of the growth are different for different sources. For a 
weak source with E = 5.8 x 10^^ erg s~^, the ionized volume at the end of the 
fast growth phase is already comparable with that at later time. However, for 
a strong source with £^ = 5.8 x 10^^ erg s~^ the ionized volume at the end of 
the fast growth phase is much smaller than that at later time. These features 
can be seen more clearly in Figure 2, which plots dlnV{t)/dlnt vs. Int for 
the same solutions in Figure 1. dlnV{t)/dlnt is the index a of the power-law 
relation y ~ t". Figure 2 shows once again the three phases of the power- law 
index. When t is small, d\nV{t)/dlD.t ~ 3, or V{t) oc t^, and the radius of the 
ionized spheres, or the I-front, ri = {3V/AnY^'^ ^ ct. Therefore, this actually is 
the fast or relativistic phase. In the transition phase, dlnV{t)/d\nt decreases 
from 3 to about 1. Finally, dl'n.V{t)/dha.t approaches to ~ 1, or ri oc this 
is the slow or non-relativistic phase. 
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Fig. 1. Ionized volume V{t) vs. time t. The source intensities are taken to be 
E = 5.8x [1039 (dash dot dot), 10^^ (dash dot), lO^^ (dash), and 10^^ (gojij ij^c)] 
erg sec-\ and N to be 1.34x [10^°, 10^^^ iq54^ ^^^^ jq56j gg^^-i^ respectively. Redshift 
is taken to be 1 + z = 10. 




1 O \ I \ \ \ \ I \ \ \ '1^1 i--^..^..L.jr'r~'i--rH-- i_ 

Int Myra 

Fig. 2. dlnV{t)/dhit vs. Int(Myrs) . The source intensities arc taken to be 
E = 5.8x [10^9 (dash dot dot), 10"^^ (dash dot), 10^^ (dash), and lO''^ (solid line)] 
erg sec~^. 
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Fig. 3. logio(ic) vs. logio-/V(sec Other parameters are the same as those in 
Figures 1 and 2. tcis in the unit of mean free flight time of ionized photons. 



2.2 E- dependence of the ionized region growth 



Wc define a time scale, t^, by d\nV{t)/d\nt\t=t^ = 2.5, which characterizes 
the transition time from the fast phase to the slow phase. The dependence 
of the time scale tc on E is plotted in Figure 3, which approximately follows 
tc oc Therefore, the transition time of the growth of the ionized volume 

is non-lincarly dependent on the emission rate of the ionizing photon E, which 
indicates that the growth of V{t) also depends non- linearly on E. 

The non-linear dependence of V{t) on E can be more prominently demon- 
strated by Figure 4, in which we plot the growth of the total ionized volume 
Vtotai(^) around one source with E = 5.8 x 10"^^ erg s~\ and those given by m 
sources with intensity 5.8 x 10*^^/m erg s~^, with m = 10^, 10^ and 10^. Here 
we assume that the ionized regions for different sources do not overlap. Figure 
4 shows that the growth rate of the ionized volume depends substantially on 
the number of sources, in spite of the fact that in all cases the total photon 
emission rate, A^, arc the same. For a single strong source of E = 5.8 x 10^^ 
erg s~^, Kotai(^) at t ^ 1 Myrs is only about 0.16, 0.10 and 0.08 of that from 
m sources with intensities E = 5.8 x lO'^^/m, erg s~^ and m = 10^, 10'^ and 
106. 

Therefore, the larger the tc, the less effective the ionization with the same rate 
A^. This is simply due to the retardation of photon propagation. In the period 
of t < tc, the I-front propagates with about the speed of light c, and therefore, 
most photons emitted later will delay their contribution to the reionization 
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Fig. 4. Evolution of Vj^otai(^)- The source intensities are taken to be £^ = 5.8 x [10^^ 
(dash dot dot), 10^^(dash dot), 10^^ (dash), and 10^^ (solid line)] erg sec^^ 



by a time tc- The longer the tc, the longer the delay. When t is larger than tc, 
the speed of the I-front is slowing down and the later emitted photons start 
to join the reionization. 

It is interesting to compare Figures 2 and 4. Figure 2 shows that the ionized 
volume growth velocity dlnV{t)/dlnt approaches to 1 at i ~ 5 Myrs for all 
sources with E < 5.8 x 10^^ erg s"-*". However, Figure 4 shows that the total 
ionized volume Kotai,ix45(i) of one source = 5.8 x 10"^^ erg s"^ at 5 Myrs 
is still much less than the total ionized volume V|;otai,mx45/m(^) of m (> 1) 
sources with E = 5.8 x lO'^^/m erg s~^. That is, the total ionized volume of 
one source £■ 5.8 x 10^^ does not catch up with the total ionized volume 
of m sources with E = 5.8 x erg s~^ even when they have about 

the same growth velocity d\nV{t)/dlnt ~ 1. This is because (iVtotai.mxsg/m/c^^ 
is larger than or equal to dVtotai,ix45 / dt for all time < t < tree, we have 

Vtotal,mx45/m(^) ~ Kotal,lx45 1-^) = /o ['^Kotal,mx45/m/'^^ ~ C^Kotal,lx45/c^^]'^^ > in 

the period t < tree, i-e. before the ionized regions approach their Stromgren 
sphere. At 1 + 2; = 10, tree — 8-6 x 10^ yrs, which is comparable with 1/H, and 

therefore part of the UV photons from strong sources may cause ionization 
when the sources have already ceased. Comparing with weak sources, the 
reionization of strong sources such as > 10"^^ erg s~^ is less effective at the 
period t < tree- 
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Fig. 5. The evolution of Vtotai(*) for the number density distribution given by eq.(2). 
The source intensities are taken to be ^ = 5.8 x [10^^ (dash dot dot), 10^^ (dash 
dot), 10"^^ (dash), and 10^^ (soUd line)] erg sec~^. 



2.3 Inhomogeneous distribution of gas 



Generally, the mass density of gas is high near the source. To study its effect, 
we assume the density distribution n{r) is given by 



n(r) ^ 



no 



+ —e 
no 



-r/R 



(2) 



where R is the size of the high density region, and Uc/uo is the density increase 
in the center of the sources. As an example, we choose {uc/uo) = 10 and R = 50 
in the unit of the mean free path of the ionization photons. The growth of 
the ionized volume from different intensity of sources is plotted in Figure 5, 
which shows the similar behavior as that in Figure 4. Quantitatively, the high 
density core of n.f/no = 10 and R = 50 plays a similar role as a sphere with 
size 500 in the unit of mean free path. We have also calculated the evolution 
using other inhomogeneous density models and have obtained results similar 
to those in Figure 5. 
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3 Clustered point sources 



3.1 Time scale with the merged ionized regions 

The UV photon sources of the first generation of stars may not be very strong, 
however, they are most likely clustered. The merging of ionized regions around 
a single source will lead to complicated configuration of the ionized regions of 
clustered sources. However, in terms of the growth rate of the ionized volume, 
the problem is simplified. The merging process of clustered ionized regions can 
approximately be decomposed into a set of two-region merging. It is similar 
to the identification of clusters by the friend-of-friend method. Therefore, we 
may reveal some common features of the merging effect on the ionized volume 
growth rate by a detailed study of the merging of two ionized regions. 

Let us consider two sources located, respectively, at {x,y,z) = (0,0, a) and 
(0,0, —a). We calculate the time-dependence of the profile of the ionized re- 
gions around the two sources. The evolution of the profiles of the ionized region 
in the p — z plane (p = ^/x'^ + y'^) is shown in Figure 6, where the variables p 
and z are dimensionless as defined in §2.1. The ionized region is still defined 
by the region in which fm{p, z, t) < 0.9. In Figure 6, the intensities of the two 
sources are taken to be ii^ = 0.5 x 5.8 x 10'^^ erg s~^, and a = 1, 10 and 100 in 
the unit of mean free path of the ionizing photons. For each case, the profiles 
are at the time t' = 10, 100, 200, 400, 600, 1000 in the unit of the mean free 
flight time. 

From Figure 6, one can see first that the configuration of the ionized regions 
is very different from the ionized sphere of a point source. In the case of a = 1, 
the ionized regions have already merged when t' < 10, and the profile of the 
merged ionized region is like a sphere. For a — 10, there are two spheres 
around the two sources when t' < 10, and they merge at t' < 100. The profile 
of the merged ionized region is no longer spherical. For a = 100, however, no 
merging occurs even at the time t' ~ 1000. This is simply because the time 
scale tc of sources with £^ = 0.5 x 5.8 x 10^^ erg s~^ is ~ 30. If the distance 
between the two sources is less than etc, the merging is realized in the fast 
phase. On the other hand, for a = 100, which is larger than etc, the merging 
time will be much larger than a/c. 

In Figure 7, we show the evolution of the ionized profiles for two sources 
located, respectively, at {x, y, z) — (0, 0, a) and (0, 0, —a) with intensities E — 
0.9 X 5.8 X 10^^ and i? = 0.1 x 5.8 x 10^^ erg s^^ Similar to Figure 6, in 
the case of a = 1, the ionized regions have merged when t < 10, and the 
profile of the merged ionized region is like a sphere. For a = 10, there are 
two spheres around the two sources when t < 10, and they merge a,t t < 100. 
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Fig. 6. The evolution of the ionized region (/hi(/0, -z, t) < 0.90) of two UV photon 
sources with intensity E = 0.5 x 5.8 x 10^^ erg s~^. The contours from small to 
large correspond, respectively, to the time t' = 10 (dash dot dot), 100 (long dash), 
200 (dot), 400 (dash dot), 600 (dash), 1000 (solid) in the unit of the mean free 
flight time, p and z are dimensionless, i.e. in the unit of mean free path of ionizing 
photons. The distance between the two sources is a = 1 (top), 10 (middle) and 100 
(bottom) also in the unit of mean free path of ionizing photons. 
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Fig. 7. The evolution of the ionized region {fmip, z,t) < 0.90) of two UV photon 
sources with intensities E = 0.9 x 5.8 x 10"^^ and 0.1 x 5.8 x 10^^ erg s~^. The contours 
from small to large sizes correspond, respectively, to the time t' = 10, 100, 200, 400, 
600, 1000 in the unit of the mean free flight time, p and z are in units of mean free 
path of ionizing photons. The distance between the two sources is a = 1 (top), 10 
(middle) and 100 (bottom) in the unit of mean free path of ionizing photons. 
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There is no merging for the case of a = 100 even when the time is as large as 
t ~ 1000. Therefore, the basic feature of Figure 7 is the same as that in Figure 
6: for two sources with the transition time tic and t2c, if their distance 2a is 
less than c[tic + he], the merging occurs quickly, while it will be very slow if 
2a > c[tic + t2c]- 

From the middle panel of Figure 7, it is interesting to see that the two ionized 
spheres have about the same size, although the intensities of the two sources 
are different by a factor of about 10. This is because in the fast phase, the 
growth of the ionized sphere radius is given by the speed of light, regardless 
of the intensity of the sources. 



3.2 Growth of the ionized regions of two sources 

For the two source case, the configuration of the ionized regions generally is 
very different from the spherical ionized region of a point source. What we 
want to show in this section is, however, that the growth of the total ionized 
volume, V{t), of two sources with intensities Ei and E2 is the same as that of 
a single source oi E — Ei + E2 if the distance between the two sources 2a is 
less than c{tic + t2c), with tic and t2c being the transition time scales of the two 
sources, respectively. If 2a is larger than c{tic + t2c) , the effect of the merging 
is small, and the growth of the total ionized volume of the two sources can 
basically be treated as two isolated ionized regions, i.e. its growth rate should 
be faster than that of the single source oi E — Ei + E2. 

Figure 8 plots the evolution of V{t) of two sources with parameters Ei and E2 
as those used in Figure 6 of the last section. The V{t) of a single source with 
intensity E — Ei + E2 is also given in Figure 8. From Figure 3, we know that 
for a = 1 and 10, 2a is less than c{tic + ^20), while for a — 100, 2a is larger 
than c{tic + t2c) ■ Figure 8 indeed shows clearly that the growth of the total 
ionized volume of a = 100 is faster than that of a = 1 and 10, although the 
source intensity of a = 100 is the same as that of a = 1 and 10. The growth 
of the ionized volume of the cases a — 1 and 10 are the same. They are also 
exactly the same as the V{t) of a single source with intensity E — Ei + E2, 
though the configuration of the ionized regions of two sources is very different 
from that of a single point source. 

In Figure 9 we plot dlnV{t)/d\nt vs. In t of the two sources in Figure 8. Similar 
to Figure 2, the evolution of dlnV{t)/dlnt generally consists of three phases: 
the fast growth phase, when t is small, dlnV{t)/d\nt ~ 3, the transition 
phase, and the slow growth phase, when t is large, dlnV{t)/d\nt approaches 
to ~ 1. For the ionized region of two sources, one can not define the I-front 
with the radius of the ionized region, because the ionized region is no longer 
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Fig. 8. Ionized volume V{t) vs. time t of two sources with the same parameters as 
those in Figure 6. a = 100, 10 and 1 are shown, respectively, by solid line, filled 
circles and crosses. The unfilled square symbols are for the single point source with 
E = El + -£/2- 




Fig. 9. d\nV{t)/d\nt vs. Int (variable t is dimensionless) of two sources with in- 
tensity E = (1/2)5.8 X 10"^^ erg s"^ at a = 1 (cross) and 100 (solid). The unfilled 
square symbols are for a single source with E = 5.8 x 10*^^ erg s~^. 

spherical. However, dlnV{t)/d\nt ~ 3 tells us that the length scale of the 
non-spherical ionized region should increase with the speed of ~ c. 

As expected, Figure 9 shows that the transition time scale tc of the case 
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Fig. 10. Ionized volume V{t) vs. time t of two sources with the same parameters 
as those in Figure 7. a = 100, 10 and 1 are shown, respectively, by solid line, filled 
circle symbols and cross symbols. The unfilled square symbols are for a single point 
source with E = Ei -\- E2. 

a = 100 is shorter than that of the cases a = 1 and 10. Figure 9 presents also 
the curve of dlnV{t)/d\nt vs. Int of the single point source with E = E1 + E2, 
which is almost identical with the curve of a = 1. Therefore, we can conclude 
that in terms of the growth of the ionized volume, two sources with distance 
a — 1 and 10 are equal to a single point source with E — Ei + E2. 

Figure 10 is similar to Figure 8, but for the two sources with parameters 
used in Figure 7. In this case, most UV photons come from source 1 with 
E^O.Qx 5.8 X 10^^ erg s"\ and the intensity of source 2 E^O.lx 5.8 x 10^^ 
erg s~^ is much smaller than that of source 1. Nevertheless, we still can see 
that the growth of the total ionized volume of a = 100 is faster than that of 
a = 1 and 10. The growth of the ionized volume of the cases a = 1 and 10 are 
the same, and it is also the same as that of the single source with intensity 
E — El -\- £'2. 

Figure 11 is similar to Figure 9, but for the two sources shown in Figure 10. 
Similar to Figure 9, the transition time scale tc of the case a = 100 is shorter 
than that of the cases a = 1 and 10. The curve of dlnV{t)/dl'n.t vs. Int for 
a = 1 is almost identical to the curve of a single point source with E — E1+E2. 

Thus, in terms of the growth of the ionized volume, once the two sources 
El and E2 have a distance < c{tic + t2c), one can replace them by a single 
source of £^ = Ei + E2. Furthermore, one may apply the merging of two 
sources to multi-sources. In other words, in a cluster of sources, we can use 
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Int Myrs 

Fig. 11. d\nV{t)/d\nt vs. Int (variable t is dimensionless) of two sources with 
intensity ^ = 0.9 x 5.8 X 10^^ and 0.1 x 5.8 x 10^^ erg s'^ at a = 1 (cross) and 100 
(solid). The unfilled square symbols are for a single source with £^ = 5.8 x 10^^ erg 
s-i. 

the so-called friend-of-friend method to identify all the sources, for which the 
distances between all the nearest neighbors i and j are smaller than c{tl + ti), 
tl being the time scale tc of source i. That is, for each point source, we plot 
a sphere around the source with radius etc corresponding to its intensity, two 
sources with overlapped spheres can be identified as a cluster consisting of the 
two sources. Repeatedly applying the method to each pair of sources, one can 
then identify a cluster consisting of all the sources of which the etc spheres are 
connected. For this cluster, the ionized volume growth can be approximately 
described by a single source ol E — Y^iEi, where Ei is the intensity of source 
i. 



4 Discussion and conclusion 

We have developed the WENO algorithm to solve radiative transfer equation 
beyond one physical dimension, which can precisely reveal the features of the 
merging of the ionized regions of UV photon sources. 

We show that the growth of the ionized volume around either one or two 
point sources generally consists of three phases: fast growth phase at the early 
evolution, slow growth phase at the later evolution, and a transition phase 
between them. The transition time tc depends significantly on the intensity 
of the photon sources E, approximately to be tc oc E^/"^. Most of the photons 
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emitted by the ionizing sources would be delayed a time tc to contribute to the 
ionization. The longer the tc the less effective the ionization. Consequently, 
linear superposition is not available in estimating the ionized volume. The 
growth of the ionized volume of multiple isolated sources Ei {i — 1,2,...) 
generally is much faster than that of a single source with intensity E — T,iEi. 

Therefore, to calculate accurately the growth of the ionized volume at t, one 
should not use eq.(l), but take into account of the retardation of UV photons. 
This effect is more substantial if the first generation of stars are highly clus- 
tered. If the ionized spheres of these sources are merging in relativistic growth 
phase, the cluster would behave as a single source with intensity equal to the 
summation of intensities of these sources. Therefore, tight clustering of UV 
sources will delay the development of reionization. 

These results may already be useful in studying the 21 cm signals from the 
reionization epoch. It has been shown that a 21 cm emission and absorption 
region will develop around a point source once the speed of the ionization front 
(I-front) is significantly lower than the speed of light (Liu et al. 2007). The 21 
cm region extends from the I-front to the front of light (r = ct); its inner part 
is the emission region and its outer part is the absorption region. Therefore, 
the 21 cm region should be formed only at the time t > tc- Since sources with 
weak intensity have small tc, while strong sources or tight clusters of weak 
sources have longer tc, sources with weak intensity produce 21 cm signal earlier, 
while strong sources, or tightly clustered weak sources produce it later. The 
results calculated from equation (1) will over-predict the ionization, and thus 
lead to less 21 cm signals. These features could be important to reconstruct 
the history of reionization with 21 cm tomography and/or cross correlation 
between redshifted 21 cm signals and emission on other bands from the center 
of the sources. 
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A Equation 



The radiative transfer equation in an expanding universe is (Bernstein, 1988; 
Qiu et al. 2006) 

dJ 1 d 

+ -V ■nJ+—{HJ)^-{k, + 3H)J + S, (A.l) 



d {ct) a duj 
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where J{t, x, u, n) is the specific intensity, a is the cosmic factor, H = a/ a, v 
is the frequency of photon, u = \nl/u. and n is a unit vector in the direction 
of the photon propagation. We take c = 1 below. The absorption coefficient 



k^, is 



(A.2) 



where the cross section a{u) — (To{vo/uY and uq — 6.3 x 10 cm . 

Assuming that the point sources i {i = l,...,n) are located at Xi, and all 
photons are emitted along the radial direction, we have 



S'(t,x,z/,n) = ^/i(t,x,z/)5(n - erj. 



(A.3) 



i=l 



where is the unit radial vector with respect to Xj. The function fi is given 

by 



fi{t,x,iy) = < 



Ei{i')/V, at X = Xi 
0, otherwise. 



(A.4) 



where Ei{i') is the energy of photons emitted from the sources i per unit time 
at frequency u. When y — > 0, fi{t,yi,v) Ei5{yii). We assume the energy 
spectrum of UV photons to be of a power law Ei^v) = Eq{vq/v)°'^ and vq is the 
ionization energy of the groimd state of hydrogen hv^ = 13.6 eV. Integration 
of E over v gives the total intensity (energy per unit time) of ionizing photons 
emitted by the source, E — E{iy)diy = Eoiyo/{a — 1). 

Since there is no photon-photon collision, the solution of eq.(A.l) can be writ- 
ten as 

n 

J = X^J,5(n-erJ (A.5) 

1=1 



and Jj satisfies following equation 
dJ 

-gj + V ■er.Ji^-k.Ji + U (A.6) 



The coupling among Jj is via the absorption coefficient ki, defined in eq.(A.2). 
The evolution of the number density of neutral hydrogen, nHi(t, x), is governed 
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by the ionization equation, 
dfm 



dt 



where the fraction of neutral hydrogen /hi = nm{t, x)/n(x), n(x) is the hydro- 
gen distribution, and Ue = n{x.) — nHi(t,x) is the number density of electrons, 
if the electrons from ionized helium can be ignored, ami is the recombination 
coefficient, FeHi is the collision ionization rate, and the photoionization rate 
r^m{t,p,z) is given by 

r,Hi(t,x) = tjdu^^^^j^aiu). (A.8) 



The evolution of the hydrogen gas temperature T, in the unit of K, is given 

by 

^ = H- n'C, (A.9) 



where U — ^nksT with be the Boltzmann constant, and 

n °° 



The relevant parameters are as follows (Theuns et al. 1998) 

1. The recombination coefficient 

ami = 6.30 x IQ-^^T-^'^T^""-^ / {1 + T"'^), (A.ll) 

where T is temperature, and T„ = T/IO". 

2. The collision ionization 

TeHi = 1.17 X io-iOrV2e-i57809.i/T(i ^ tI''')-\ (A.12) 

3. The cooling. Since only the recombination cooling is important, we have 

C = 8.70 X 10-'^Ti/2y3-°-2(l + r°-^)-i[l - fnif (A.13) 
+ 1.42X 10-2^tV2[i_/hj]2 
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+ 2.45 X 10-21rl/2e-157809.1/T(i ^ j.V2)-l(i _ f^^^f^^ 

+ 7.5 X 10-^^e-^^«^^«/^(l + T^^Y\1 - /hi)/hi 

where Tn = T/10". The terms on the r.h.s. of eq.(A.13) are, respectively, 
the recombination cooling, the coUisional ionization cooling, the coUisional 
excitation cooling and bremsstrahlung. Both H and C are in the unit of ergs 
cm^ s~^. 



B Two sources 



Consider the case of two sources. Using cyhndrical coordinate {p,6,z), the 
two sources are assumed to be located at x = x± = (0,0, a), with e,.^ = 



(p/ Y + (-2 + o-Yi 0, /yp^ + (2; + a)^). The corresponding radiative trans- 
fer equations eq.(A.6) are 



9t pdp \^yy + (7+«0' / \^Jp^ + {z^aY 
= -KJ± + /±, 

where J±{t, p, v) are the specific intensities for sources X-|-. 

Instead of adding the source term f± in the r.h.s of eq.(B.l), equivalently, we 
impose a boundary condition 

lim 47rr|j±(t,p,^,i/) = E^{v), (B.2) 

p— »U,2— +±a 



where r± — p"^ -\- {z^ of, — Eqj^{vq/v)°' is the energy of photons 

emitted from the sources per unit time at frequency u. 

The absorption coefficient in eq.(B.l) ky is defined in eq.(A.2), with nHi(t, p, z) e 
^/hi(^, Pi z) governed by the ionization equation eq.(A.7). Here the photoion- 
ization rate r^Hi(^,P, -2) in eq.(A.7) is given by 



r^ui{t:P:Z) = J dv — a[v). (B.3) 



The kinetic temperature of the baryon gas is determined by eq.(A.9) with 

00 

H^nm j{J+{t,p,z,u) + J_{t,p,z,u))a{u)'^—p-diy. (B.4) 
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C The numerical algorithm 



In this paper, we are solving the system of equations (B.l), (A. 7) and (A. 9). 
To apply the WENO algorithm, we rewrite eq.(B.l) into a conservative form 
as 

dJ± d , p ^ , d ,z ^ a ^ , 1^ ,^ /^-.N 

+ ^[—J±] + ^^^^± = J± - KJ±, C.l 

at op r± oz r± r± 

where r± — + (2; ^ a)^. We use the WENO algorithm to approximate the 
spatial derivatives in eq.(C.l) (Qiu et al. 2006). 

In the numerical implementation, it is convenient to introduce the dimension- 
less variables f, p', z', a', u', J' by rescaling t' — cnaot, p' — naop, z' = naoz, 
a' — naoa, v' — vjvQ^ and J'^dv' — ■j^^J±di'. Therefore, t' and r' are respec- 
tively, the time and distance in the units of mean free flight time and mean 
free path of ionizing photon huQ in the non-ionized background hydrogen gas 
n. For the ACDM model, n = 1.88 x 10~^(1-|-^)^ cm~^, where z is the redshift, 
t = 0.89(1 + z)-H' Myrs and r = 0.27(1 + z)-^r' Mpc. Then the system of 
equations (C.l), (A. 7) and (A. 9) can be rewritten as the following system 

9J± d p' d z'Ta' w 1 1 1 . .X 

c^o^ = o^Hiiil - fmf - ^fm - TeHiil - Mfm (C.3) 
dt n 

-caokB-^^H-C (C.4) 



with ^^^^ given by 



-r,Hi{t,p,z)^ J^^^du', (C.5) 



H given by 



H = huofm [ (j; + J'S-^dv', (C.6) 



1 



and anii-i ^eHi and C given by equations (A. 11), (A. 12) and (A. 13) respec- 
tively. 



20 



To solve the radiative transfer equation (C.2), wc adopt the fifth-order finite 
difference WENO scheme, which was designed in (Jiang & Shu 1996), coupled 
with the third order TVD Runge-Kutta time discretization for the system of 
equations (C.2), (C.3) and (C.4). The multi-time-scale strategy (Qiu et al. 
2007) and the adaptive time step strategy (Liu et al. 2007) are used to save 
the increased computational cost introduced by the stiffness of the equations 
(C.3) and (C.4). The numerical algorithms are implemented as describe below. 
For the sake of simplicity, we drop the prime in the notations. For example, J 
means J' hereafter. 

• The computational domain and computational mesh: 

The computational domain is (p, z, v) G [0, Pmax] x [-z„iax, z„,ax] x [1, i^rnax], 
where Pmax and z^ax are chosen such that J{t,p,z,p) ~ for p > Pmax, 
\z\ > Zmax or 1/ > Vmax- In our Computation, pmax is taken to be greater than 
the final computational time t, Zmax is taken to be Pmox + a and Vmax — 10^. 

To avoid tj^ in cq. (C.2) to become 0, we design the computational mesh, 
such that x-t arc located at the center between grid points. The mesh sizes in 
the p- and z- directions are set to be the same, i.e. Ap = Az, to preserve the 
spherical symmetry property of the numerical solution before the merging 
of the two sources. The mesh in the v direction is taken to be smooth but 
not uniform. Specifically, the mesh sizes are designed as the following 

= " Ap = A^; A{ = logs ^max/N,,] 

with A^^^ being the number of mesh points in [0, a] in the 2;-direction, and 
being the number of mesh points in the i/-direction. The computational 
mesh is 

p,= (^i--)Ap, i = l,2,...,Ar^, 

Zj=iAz, j = 1,2, A^2, 
z/fe = 2^*, with = kA^, k = 1, 2, N^. 

• The WENO method in approximating the spatial derivatives: 

The approximation to the point values of the solution J±{t"', pi, Zj, i/^), de- 
noted by J±^ij^k, is obtained with a dimension by dimension approximation 
to the spatial derivatives using the fifth order WENO scheme (Jiang & Shu 
1996). Taking ^(^^J+) as an example, the approximation is performed 
along the 2;-line with fixed p, and i/^: 

where the numerical flux hj_^i is obtained with the following procedure. 
When the "wind direction", namely the coefficient ^^"'"^^+^ — a is positive 
at the mesh boundary, we use a left-biased stencil in reconstructing the 
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numerical flux h,,i as described in detail below. When the "wind direction" 
is negative, we use a right-biased stencil to obtain the numerical flux hj^i, 
following a mirror symmetry reconstruction with respect to j + | as that of 
the left-biased stencil. When the coefficient ^■'^^■'^^ — a = 0, which actually 
will happen due to the way we design our mesh, the numerical flux is simply 
set to be 0. 
We denote 

hj = J(r, Pi, Zj, i/fe), j = -2, -1, + 2 

where n, i and k are fixed. The numerical flux from the regular WENO 
procedure is obtained by 



ere ^^+1/2 three thii 

by 



where ^^+1/2 three third order fluxes on three different stencils given 







7, 

- 6^.-1 


6 


^^^■+1/2 - 


6 


5^ 

-1 + -h, 


+ 3^j+i 


^^^■+1/2 - 


1^ 

3/^.+ 


5, 





and the nonlinear weights ujm are given by 

li 



COm — ~i 



with the linear weights 7; given by 

_ 1 _ 3 _ 3 

71 - 72 - 73 - 

and the smoothness indicators (5i given by 

A = ^ (/i,-2 - 2/i,_i + /i,)' + J (/i,_2 - 4/i,_i + 3/i,)' 

13 1 

" 12 ^^^"^ ~ ^ ^^'^^''^ ^ 4 ^^^"^ ~ 
13 1 

" 12 ~ ^ ^^^^^^ ^ 4 ~ ^ ^^^^^^ ■ 

£ is a parameter to avoid the denominator to become and is taken as 
e = 10~^ times the maximum magnitude of the initial condition J in the 
computation. 
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• Time integration: 

To evolve in time, we use the third order TVD Runge-Kutta method 
(Shu & Osher 1988). For systems of ODEs Ut = L{u), the third order 
Runge-Kutta method is 

+ AiL(ii",r) (C.8) 
1,(2) = ^^" + 1(^^(1) +AiL(^/«)) (C.9) 

3 3 

The difficulty of the direct implementation of the Runge-Kutta method 
hes in the stiffness of equations (C.3) and (C.4). Especially for the strong 
source, one needs a very small time step At, as small as 10~^, to guarantee 
the stability of the numerical scheme, therefore the computational cost for 
long time integration is huge. In this paper, we adapt the multi-time-scale 
strategy (Qiu et al. 2007) and the adaptive-time-step strategy (Liu et al. 
2007) to evolve the system of equations in time. We refer to (Qiu et al. 
2007) and (Liu et al. 2007) for the details of its implementation. 

• Numerical boundary condition: 

• In the p-direction, 
at p = 0, 

J±,-i,j,k — J±,i+l,j,k, ^ = 0, 1, 2, 

^ P — Pmax-i 

J±,Np+i,j,k = 0, i = 0, 1, 2. 

• In the z-direction, 

J±,i,T(.N^+i),k = 0, z = 0, 1, 2 

• Around the point sources x-t, according to eq.(B.2), 

1 E 

with r±^ij — ^ Pi + {zj ^ a)2. is a small number depending on the mesh 
size. Tg is bigger when the mesh is coarser. 

• Parallel computing: 

The computational cost in solving the system of equations (C.2), (C.3) 
and (C.4) is quite large for this 3-dimensional (2-D in the physical space and 
1-D in the frequency space) time dependent problem. Parallel computing 
with mpif77 is used to speed up the computation. 
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